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The probability distributions of the order parameter for two models in the directed percolation uni- 
versality class were evaluated. Monte Carlo simulations have been performed for the one-dimensional 
generalized contact process and the Domany-Kinzel cellular automaton. In both cases, the density 
of active sites was chosen as the order parameter. The criticality of those models was obtained by 
solely using the corresponding probability distribution function. It has been shown that the present 
method, which has been successfully employed in treating equilibrium systems, is indeed also useful 
in the study of nonequilibrium phase transitions. 
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I. INTRODUCTION 

The theoretical treatment of nonequilibrium dynamic 
systems has been given a great deal of attention dur- 
ing the last decades. However, even the simplest models 
have not yet reached the same level of understanding as 
their equilibrium counterparts. Excellent resources on 
this subject can be found in Refs. [l|-[3| and references 
therein. 

One of the most important models in nonequilibrium 
physics is directed percolation (DP), that is widely con- 
sidered as an analogue of the Ising model for nonequi- 
librium phase transitions. DP has been used to simulate 
a large variety of problems, including flow of a liquid 
through a porous medium, electric current in a diluted 
diode network, and reaction-diffusion processes. Another 
interesting model is the so-called contact process (CP), 
which was first introduced by Harris as a nonequilib- 
rium toy model to study epidemic spreading. In the stan- 
dard CP, each site of a lattice can be active, representing 
an infected individual, or inactive, corresponding to a 
healthy person. The system evolves in such a way that 
only one site is updated at a time. In a d-dimensional 
hyper-cubic lattice, the annihilation rate /i means that 
an active (occupied) site becomes inactive (vacant) at 
rate /z, independent of its neighbors. A vacant site turns 
to occupied at a creation rate proportional to the frac- 
tion of occupied neighbors, n/2d, where n is the num- 
ber of occupied nearest neighbors. Thus, an inactive site 
surrounded only by inactive neighbors remains inactive. 
Once all sites are vacant, the system becomes trapped in 
that state, which is known as frozen or absorbing state. It 
is well known that the CP undergoes a second order phase 
transition from an active to a frozen (absorbing) phase at 
some critical /i c . For annihilation rates /x > /i c , the only 
quasi-stationary state is the absorbing one, while for suffi- 
ciently small /i a finite fraction of sites remains active. In 
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spite of its simplicity, no exact results for fi c are known, 
even for one-dimensional lattices. Mean field approxima- 
tion, Monte Carlo simulations, and series expansion are 
the most used techniques JEHU. The best estimate for the 
critical point in one dimension 0, Q is /i c = 0.303228(2). 

A generalized version of the CP is obtained by consid- 
ering different creation rates. For instance, in one dimen- 
sion, let us call £ the creation rate at a given site with 
exactly one occupied neighbor. The creation rate with 
two occupied neighbors is set to 1, while annihilation oc- 
curs at rate zx. Standard CP corresponds to £ — 0.5, 
while C = 1-0 is known as the A-model[9(. This family 
of processes was proposed by Durret and GriffeathflG] 
and has received some attention later on0, [ll|. Anal- 
ogous to the usual CP, this generalized contact process 
(GCP) also shows a continuous phase transition from the 
active state to the absorbing state for infinite systems. A 
different generalization of the contact process, although 
not discussed in the present study, considers more than 
one absorbing state [12J and is also a subject of current 
interest [lj, [lj| . 

Another irreversible model that describes a nonequi- 
librium phase transition from an active to an absorbing 
state is the probabilistic Domany-Kinzel cellular automa- 
ton (DKCA)[l5j]. It is represented in a one-dimensional 
lattice containing N sites that can be either empty or 
occupied. The state of each site i at time t + 1 depends 
only upon the state of the two nearest neighbors i — 1 and 
i + 1 at time t. In contrast to the GCP, in the DKCA all 
sites are updated simultaneously. Denoting the occupa- 
tion variable by cr,; , one can define the conditional proba- 
bilities P[(Ti(t + l)\ai-i(t),<T i+ i(t)], where <7j(i) =1(0) if 
the site i is occupied (empty) at time t. Besides the ac- 
tive/absorbing transition, this model exhibits a damage 
spreading transition in the active phase, from a chaotic 
to a non-chaotic region [iH 

Regarding universality classes, it is believed that most 
models having absorbing state transitions belong to the 
directed percolation universality class since some essen- 
tial features like short-range interactions and transla- 
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tional invariance are fulfilled [ 1 8M2 Oj . Thus, the absorbing 
transitions in the GCP and the DKCA models are in the 
DP universality class From an experimentalist point 
of view, this class was recently observed in turbulent liq- 
uid crystals |22l l23f. 

In what concerns physical thermodynamic quanti- 
ties, the order parameter distribution function has been 
proven to be an important tool for stud ying a large va- 
riety of subjects, as magnetic systems |24l428j . the liquid- 
gas critical point (29|, the critical point in the unified the- 
ory of weak and electromagnetic interactions I30l , and the 
critical point in quantum chromodvnamics [3 lj . To the 
best of our knowledge, even with all those applications, 
the method of the order parameter analysis did not re- 
ceive much attention in nonequilibrium physics. In the 
present work, we explore the probability distribution of 
the order parameter of the GCP and the DKCA in or- 
der to analyze the DP universality class. Results show 
that the method is indeed quite reliable to study nonequi- 
librium phase transitions, and we hope that it could be 
generalized and applied to other dynamic models. 



II. APPROACH 

For the specific cases of the GCP and the DKCA, the 
order parameter can be chosen as the density of active 
sites, namely p — Y2i=i °i) where N is the total num- 
ber of sites and the occupation variable cr.; is equal to 1 
(0) if site is active (inactive). In the infinite-size limit, p 
vanishes in the absorbing phase. In finite-size systems, 
the density p is a fluctuating quantity, characterized by 
the probability distribution P(pJ . Analogous to the usual 
finite-size scaling assumptions [32j, one then expects that, 
for a large finite system of linear dimension L at the crit- 
ical point, P(p) takes the form 



P(p)=bP*(p), 



(1) 



where b = b( S L^/ v - L , j3 and i/± are the critical expo- 
nents of the density of active sites and the correlation 
length, respectively, p = bp, &o is a non-universal con- 
stant, and P*(p) is a universal scaling function. For 
the DP universality class, one has /? — 0.276486(8) and 
z/_L = 1.096854(4) [33j. Scaling functions, such as that 
given by Eq. ([I]), are characteristic of the corresponding 
universality class. Systems belonging to the same univer- 
sality class share the same P* scaling function and thus, 
from the precise knowledge of P* (p) , one can characterize 
critical points and also identify universality classes. 

The most efficient way to compute the probability dis- 
tribution P(p) is probably through Monte Carlo simula- 
tions. In equilibrium systems, P{p) corresponds to the 
fraction of the total number of realizations in which the 
order parameter reaches the specific value p. In absorb- 
ing state systems, obtaining that distribution is a more 
complicated issue, since the active stationary distribution 
only appears in the infinite-size limit. For finite lattices, 



as the system always becomes trapped in the absorb- 
ing state, one can only evaluate the quasistationary (QS) 
distribution 

Let us briefly review the definition of the QS distri- 
bution. By denoting n as the number of active sites 
(n = corresponds to the absorbing state) and P n (t) 
as the probability of having exactly n occupied sites 
at time t, the survival probability can be obtained by 
Ps{t) = J2n=i p n(t) = 1 - Po{t). As t -»• oo, it is ex- 
pected that P n , normalized by the survival probability 
p s (t), remains time-independent]!]. A procedure to com- 
pute the QS distribution is to restrict averages over the 
surviving realizations only, i. e., after performing a large 
sample of independent realizations, the average value of 
some physical quantity at time t is taken over the real- 
izations that did not reach the absorbing state at that 
time. At long times, as the number of surviving samples 
decays, this mechanism suffers from large fluctuations. 

A more effective way to compute the QS distribution 
was proposed by Tome and de 01iveira[34j. It consists 
in creating a particle in the finite system whenever the 
absorbing state is going to be reached. This procedure 
is equivalent to forbid the last particle to be annihilated 
and thus the density of active sites p is always non-zero. 
In the thermodynamic limit this perturbation was found 
to be irrelevant, as shown in Ref . [33| . The same authors 
have also proposed a conserved contact process in which 
p is constant and the absorbing state is eliminated [35j. 
The model can be seen as the CP version in an ensemble 
of fixed particle number and its properties, in the ther- 
modynamic limit, are identical to those of the ordinary 
CP. The equivalence between both ensembles was shown 
by Hilhorst and van Wijland 36]. 

Another powerful method to obtain QS distributions 
was proposed by de Oliveira and Dickman(37j. It con- 
sists in storing a list with M non-absorbing configura- 
tions that the system has visited previously (typically 
M ~ 10 3 - 10 4 ). The list is updated with probability 
p r (usually p r ~ 10~ 3 — 10~ 2 ), which means that a con- 
figuration from the list is replaced by the current con- 
figuration with probability p r . During the simulation, if 
an absorbing configuration is imminent, it is replaced by 
another one, randomly chosen from the list. This pro- 
cedure is used in the present work, with M = 2000 and 
p r = 10~ 3 in most cases. 

Regarding the simulation details, we have simu- 
lated the generalized contact process (GCP) in one- 
dimensional lattices with periodic boundary conditions 
and sizes L varying from 80 to 640, and up to 3200 in a 
few cases. Different starting configurations were tested, 
with an initial density of active sites pi varying from 
0.2 to 1, and the QS distribution was found to be in- 
dependent of pi, within the error bars. For each lattice 
size, simulations of 10 — 1000 samples with 10 6 — 10 7 
Monte Carlo steps per sample were performed. Transi- 
tion rates are schematically represented in Table Q] Ac- 
cording to those transition rates, the time evolution can 
be described as following [38j. 
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• Choose a site i randomly. 

• Choose a process (creation or annihilation) : 

— for C < 1 ; choose creation with probability 
j^j2 and annihilation with probability 

— for £ > 1: choose creation with probability 

and annihilation with probability -^-^ ■ 

• If site i is vacant (&i = 0) and creation was chosen, 
one should define n = Oi-\ + <Ji+\. Again, we need 
to consider both situations: 

— for £ < 1: creation occurs with probabilities 0, 

and 1 for n equal to 0, 1, and 2, respectively; 

— for £ > 1: creation occurs with probabilities 
0, 1, and 1/C for n equal to 0, 1, and 2, re- 
spectively. 

• After choosing N sites, increase time by one unit. 



The same behavior was observed for all other values of £ 
as well as for the DKCA. Figure [2] shows the normalized 
probability distribution function of the DKCA for p 2 = 
and two values of p\ . 




P 



TABLE I. Transition rates in the GCP. Open (filled) circles 



indicate vacant (occupied) sites. 

From To Rate 

o o o o • o 

• o o • • o C, 
o o • o • • £ 

• o • • • • 1 
• O fl 



FIG. 1. Probability distribution of the order parameter in the 
GCP for C = 0.5. (a) /i = 0.295 (/i < /i c ) and (b) fi = 0.305 
(/i > He)- Error bars were omitted for better visualization 
and are typically around 0.3%. 



The one-dimensional DKCA was simulated on lattices 
with up to 3200 sites and averages were done over 10 3 
samples with 10 6 Monte Carlo steps per sample. The 
transition probabilities P[<7j(t + l)|(jj_ 1 (t), <j i+1 (t)} were 
P[1|0,0] = 0, P[1|0,1] = P[1|1,0] = pi, and P[l|l,l] = 
p 2 - Naturally, P{0\a i _ 1 ,a i+1 } = 1 - P[l\a i _ 1 ,a i+ i}. As 
already mentioned, in the DKCA all sites are updated 
simultaneously. 



III. RESULTS 
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Following the mechanism proposed by Martins and 
Plascak[2lJ, we analyzed the function P*{p) to get an es- 
timate of the critical point. As depicted in Fig.[T]one can 
see that, as the lattice size increases, the peak of the func- 
tion moves to the right for n = 0.295, and it goes to the 
left for [i = 0.305. From a different point of view, let us 
consider the function P*(p) for L = 400 and \i\ = 0.295 
as shown in Fig. QJi. The same distribution shall be ob- 
tained for a larger lattice (say for instance L = 800) at 
a different rate /i2 in such a way that /i2 > Mi- O n the 
other hand, if we consider as reference the distribution 
for L = 400 and /ii = 0.305, we will have the same dis- 
tribution for a larger lattice at /i 2 < fJL\. This suggests 
that the critical \i c is in the range 0.295 < /i c < 0.305. 



FIG. 2. Probability distribution of the order parameter in the 
DKCA for p 2 = 0. (a) pi = 0.809 and (b) pi = 0.810. Error 
bars were omitted for better visualization and are typically 
around 0.3%. 

In order to obtain a better estimate of the criti- 
cal rate /x e for the infinite lattice, one can proceed as 
By using a reference distribution func- 
C, and /i, one can vary /i for a 



following [27 1 . 
tion for a given L, 
different lattice size until a distribution that collapses 
into the reference one is obtained. For instance, Fig- 
ure^ shows the normalized probability distribution for 
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L = 640 at C = 0.3 and p, = 0.19080(3), considered as 
reference. For L = 320, the same distribution is ob- 
tained at p = 0.19070(5). For L — 160, the correspond- 
ing value of fi is 0.19050(5), while for L = 80 one has 
p = 0.19020(10). All these four distributions are depicted 
in Fig. (3^. Each one of those values of p gives an estimate 
for the pseudo-critical p^ for that lattice size. Since one 
expects that the difference \pl — Mc| scales as L~ 1 / v± , 
where i>±_ is the correlation length critical exponent, a 
finite-size scaling analysis can be performed to estimate 
the critical values of the infinite system. In Fig. (SJd, 
one has a plot of p L vs. L^ 1 ^ 1 - , with vy = 1.096854(4) 
(Ref. [33]). Each row in Table [TT1 contains the values of 
\xl that lead to the same distribution function for each 
£ as well as the extrapolated value of p c , obtained from 
the finite-size scaling technique. If another distribution 
is used as reference, a different set of pl is obtained (as 
shown in Table HlT| . providing another estimate for p c . 
A similar analysis for the DKCA with p2 = is also 
depicted in Fig. [31 with the corresponding data repre- 
sented in Tab. IIVI A finite-size scaling analysis leads to 
a critical value p\. c — 0.80932(1), which has even higher 
precision than previous works (pi, c = 0.811(1) from [39j, 
and p 1<c = 0.8095(3) from [H). 



2 4 6 0.005 0.01 0.015 0.02 




FIG. 3. Graphs obtained from the data presented in Tabs. HT1 
Mil and IIVI Error bars were omitted for better visualization, 
(a) Normalized probability distribution function for the GCP 
with ( = 0.3. The plot is, in fact, a superposition of four 
curves, (b) Corresponding finite-size scaling analyses for the 
GCP. (c) Normalized probability distribution function for the 
DKCA with p2 = 0. The plot is, in fact, a superposition of 
five curves, (d) Corresponding finite-size scaling analyses for 
the DKCA. 

An often used technique to study the criticality of the 
DP universality class consists in evaluating the moment 
ratio m = (p )/{p) 2 ^^- This quantity is analogous to 
the reduced fourth cumulantflij and reaches a universal 
value at the critical point. Thus, the curves for m(p 1 L) 
cross near p c for different L. Figure @] illustrates m(p, L) 
for different values of £. To compare the results obtained 



TABLE II. Each row contains the annihilation rates /il at 
which the distribution functions are the same. Errors in 
parentheses affect the last digits. For £ = 0.01 the lattice size 
L = 240 was used instead of L = 80, giving p L = 0.00834(1). 
The last column shows the extrapolated values for /x c (C)- 
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L = 640 


L = 320 


L = 160 


L = 80 


L — > oo 


u . U 1 










u.uuouu^o ) 


0.02 


0.01638(1) 


0.01630(1) 


0.01617(2) 


0.01605(5) 


0.01641(3) 


0.05 


0.03790(5) 


0.0377(1) 


0.0374(1) 


0.0369(1) 


0.03804(3) 


0.1 


0.0709(1) 


0.0706(1) 


0.0701(1) 


0.0694(2) 


0.07111(7) 


0.2 


0.1320(1) 


0.1314(1) 


0.1305(1) 


0.1289(2) 


0.13247(7) 


0.3 


0.19080(3) 


0.19070(5) 


0.19050(5) 


0.1902(1) 


0.19090(2) 


0.4 


0.24740(2) 


0.24720(3) 


0.24685(5) 


0.2465(1) 


0.24750(8) 


0.5 


0.3021(1) 


0.3010(2) 


0.2992(2) 


0.2958(3) 


0.3031(1) 


0.6 


0.35660(5) 


0.3552(1) 


0.3528(3) 


0.3488(5) 


0.3578(2) 


1.0 


0.5735(1) 


0.5728(2) 


0.5718(3) 


0.5700(5) 


0.5740(1) 


2.0 


1.1030(1) 


1.1023(1) 


1.1015(3) 


1.1001(8) 


1.1034(1) 



TABLE III. Sets of /il obtained from another reference dis- 
tribution. Again, for £ = 0.01, L = 240 was used instead of 
L = 80, giving y. L = 0.00876(1). 



c 


L = 640 


L = 320 


L = 160 


L = 80 


L — ► oo 


0.01 


0.00870(1) 


0.00873(1) 


0.00882(2) 




0.00865(3) 


0.02 


0.01655(1) 


0.01659(2) 


0.01672(4) 


0.01705(10) 


0.01643(3) 


0.05 


0.03820(1) 


0.03827(3) 


0.03840(6) 


0.03870(10) 


0.03810(3) 


0.1 


0.07150(2) 


0.07175(6) 


0.0722(1) 


0.0730(1) 


0.07125(8) 


0.2 


0.13300(5) 


0.13330(5) 


0.1339(1) 


0.1350(2) 


0.13264(7) 


0.3 


0.19110(5) 


0.1913(1) 


0.1915(1) 


0.1921(1) 


0.19094(4) 


0.4 


0.24780(5) 


0.24800(5) 


0.2483(1) 


0.2492(2) 


0.24752(7) 


0.5 


0.3040(1) 


0.3046(3) 


0.3058(4) 


0.3080(5) 


0.3033(1) 


0.6 


0.3585(1) 


0.3588(1) 


0.3592(1) 


0.3607(4) 


0.3580(2) 


1.0 


0.5750(1) 


0.5757(1) 


0.5770(3) 


0.5800(5) 


0.5740(1) 


2.0 


1.1042(1) 


1.1045(1) 


1.1053(3) 


1.1072(8) 


1.1035(1) 



by using the probability distribution function to those 
coming from the crossings of the moment ratio, Table 
IVl shows the critical p, c achieved by both methods. The 
results that come from the probability distribution are 
the mean value of the extrapolated p c depicted in Tables 
Inland [HI! for the generalized contact process and in Table 
IIVI for the Domany-Kinzel cellular automaton. 



IV. CONCLUSIONS 

In focusing on the study of the probability distribution 
of the order parameter in systems that do not obey the 
detailed balance, this work considered two different mod- 
els in the directed percolation universality class. The gen- 
eralized contact process and the Domany-Kinzel cellu- 
lar automaton were investigated. The criticality of both 
models was obtained by using the probability distribution 
of the order parameter itself and the results showed that 
this approach is also powerful to study nonequilibrium 
phase transitions, regarding their universal and nonuni- 
versal aspects. In general, the critical values obtained 
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TABLE IV. Results for the DKCA with P2 = 0. Each row shows the values of pi t L at which the normalized probability 
distribution is the same. These data were used to plot Fig. [5Ji and give the estimate of the critical pi jC in the last column. 
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TABLE V. Comparison between the critical values from two 
different methods. 
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c 


/i c from 
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1.5 




cumulant crossings 


probability distributions 




1.4 
1.3 


0.01 


0.00865(3) 


0.00863(3) 




1.2 


0.02 


0.01648(3) 


0.01642(2) 




1.1 


0.05 


0.03810(6) 


0.03807(4) 


r 




0.1 


0.0713(1) 


0.07118(8) 




1.4 


0.2 


0.1326(1) 


0.13256(8) 




1.3 


0.3 


0.1909(1) 


0.19092(4) 






0.4 


0.2476(1) 


0.24751(8) 




1.2 


0.5 


0.3032(1) 


0.3032(1) 




1.1 


0.6 


0.3582(2) 


0.3579(2) 






1.0 


0.5742(3) 


0.5740(1) 






2.0 


1.1038(5) 


1.10345(15) 


0.6 






DKCA 


0.8093(1) 


0.80932(1) 



FIG. 4. Moment ratio m (as denned in text) as a function 
of the annihilation rate /x. Error bars when not shown are 
smaller than the symbols. 



from the present method have higher precision than the 
values from the crossings of the moment ratio. In addi- 
tion, this work has provided an accurate estimate for the 
critical point in the Domany-Kinzel cellular automaton 
with p2 = 0. To the best of our knowledge, the present 
approach, using just the probability distribution of the 
order parameter as expressed in Eq. (JlJ, was applied 
to nonequilibrium systems for the first time. We believe 
that these results will spread the treatment of other dy- 
namic systems within the present approach. Applications 
to damage spreading transitions, that are supposed to be 
in the same universality class, are now in progress. 
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